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I. INTRODUCTION 

In the functional formalism^ the general quantum 
mechanical amplitude A(a,b;T) = (b\e~ TH \a) is given 
in terms of a path integral which is simply the N — > oo 
limit of the expression 

A N (a, b; T) = (7^-) ' J % ' • • d QN-i e~ S " . (1) 

The euclidean time interval [0, T] has been subdivided 
into N equal time steps of length — T/N, with qo = a 
and qN — b. Sn is the naively discretized action of the 
theory. In this paper we will look at theories with action 
of the form 

s =£ dt Q^+^fo)) • ( 2 ) 

Note that we use units in which h and particle mass have 
been set to unity. The naively discretized action is in this 
case simply 

N-l t a \ 
Sn = E {^+^V(q n ) \ , (3) 
n— ^ ' 

where 5 n = q n+1 - q n , and q n = \{q n+ x + q n ). Key inves- 
tigations regarding numerical evaluation of path integrals 
were presented in the reviews of Barker and Henderson 3 , 
Kalos and Whitlock 4 , and Ceperley^, as well as in the pa- 
pers by Pollock and Ceperley^ and Barker 7 . A modern 
and extensive reference on the subject of path integrals is 
given in the latest edition of the textbook by Kleinert 8 . 

As we can see, the very definition of path integrals 
makes it necessary to make the transition from the con- 
tinuum to the discretized theory. This discretization, 
however, is far from unique. In fact, the details of the 
discretization procedure are extremely important both 
for analytical and numerical treatment of path integrals. 
This dependence on discretization procedure has been 
one of the principle impediments to creating a consistent 



mathematical theory of path integration. On the numeri- 
cal side, this manifests itself in the fact that path integral 
simulations remain notoriously demanding of computing 
time - so much so that certain path integral calculations 
serve as benchmarks for new generations of supercom- 
puters. 

For this reason let us note that we have the freedom to 
make two important choices that will not affect the final 
result, i.e. the continuum amplitude we seek to calculate. 
First, we have the freedom to choose the point in the in- 
terval [q n ,q n +i] in which to evaluate the potential V. 
It is well known that different points correspond to dif- 
ferent ordering prescriptions in the operator formalism. 
The choice of the middle point q n is the most common 
one. It corresponds to the symmetric or Weyl ordering of 
operators p and q, and so always leads to a hcrmitcan ex- 
pression for the hamiltonian H . Two other prescriptions 
are also often used. The left ordering prescription evalu- 
ates the potential at q n , the left boundary of the above 
interval (in the operator formalism this corresponds to 
taking the p's to the left of the g's in all the products 
that appear in the hamiltonian). Similarly, one defines 
the right ordering prescription. Although they lead to 
somewhat simpler looking expressions, the left and right 
prescriptions do not in general give hermitean hamilto- 
nians. Let us note, however, that the class of theories 
given by eq. J5J has a hamiltonian that is the sum of 
a p-dependent kinetic term and a g-dependant potential 
term and so has no ordering ambiguities. In this case 
different prescriptions lead to the same continuum am- 
plitude - the discretized amplitudes do differ, but they 
tend to the same continuum limit. 

The second, and more important freedom related to 
our choice of discretized action has to do with the free- 
dom to introduce additional terms that explicitly vanish 
in the continuum limit. We will designate such discrete 
actions as effective actions. For example, the term 

N-l 

J2 e N6lg(q n ), (4) 

n=0 



2 



where g is regular when ejv — > 0, does not change the 
continuum physics since it goes over into e 2 N dt q 2 g(q), 
i.e. it vanishes as e N . Although such additional terms 
do not change the continuum physics they do affect the 
speed of convergence to that continuum limit. 

The aim of this paper is to give a detailed exposition 
of a systematic analysis that leads to the best solution 
in the class of all equivalent effective actions, e.g. the 
effective action that leads to the fastest convergence of 
the associated discrete amplitude An to the continuum 
expression. The uncovered speedup in the path integral 
algorithm is a direct consequence of new analytical input 
that has come from the study of the relation between 
discretizations of differing coarseness of the same theory. 
We have given a brief presentation of these ideas in a 
recent papeA 

The calculations that will be presented turn out to be 
simplest in the mid-point prescription. Before we proceed 
with them it is useful to spend the next section in a 
brief overview of known results dealing with the speed of 
convergence to the continuum limit. 

II. BRIEF OVERVIEW OF KNOWN RESULTS 



Although not related to the central investigation in this 
paper, let us present a proof of eq. I© as an illustration 
of the use of the Trotter formula 

e A+6 = lim (e A ' N e & l N ) N . (10) 

Using it we can easily show the validity of the formal 
expression A-n(cl, b; T) le ^ = A N (a,b;T)"9 ht . On the 
other hand, from eq. (Q) we find that A N (a,b;T)"9 ht = 

e tN(V(a)-V(b)) AN ^ b . T yeft_ Ag a result we gee that the 

left hand side of eq. © is simply the average of An and 
A-n and so has an expansion in even powers of I/TV. 

We end this section by presenting the result obtained 
by Takahashi and Imadaii and independently by Li and 
Broughton 12 . In these papers the authors used a general- 
ized formic of the Trotter formula to increase the speed of 
convergence of the discretized partition function. Their 
final result is a derivation of a formula for the effective 
potential V eff — V + ^e^iV') 2 . The authors showed 
that by using this effective potential (in the left prescrip- 
tion) one gets 

Z N (P) eff = Z(/3) + 0(l/N 4 ) . (11) 



In this section we will compare the speed of conver- 
gence of several different prescriptions to the continuum 
limit. The naive mid-point prescription satisfies 



A N {a,b-TY 



A(a,b;T) + 0(l/N) , 



(5) 



for all a and b. By naive we mean that we use the naively 
discretized action given in eq. J2J). On the other hand, 
in the naive left prescription the amplitude for a — > a 
converges much faster 



A N (a, a; T) left = A(a, a; T) + 0(l/N 2 ) 



(G) 



This behavior can be easily shown both analytically and 
numerically. We note in passing that this is strongly 
related to the well known result for the partition func- 
tion evaluated using naive discretization in the left 
prescription^, which follows directly from the above am- 
plitude by integrating over a (to get the trace) and writ- 
ing the time of propagation T as the inverse temperature 

a, 



Z N ((3) left = Z(f]) + 0(l/N 2 ) 



(7) 



However, going back to the language of amplitudes, it is 
also easy to show that the amplitudes for different initial 
and final states converge slower, i.e. for a ^ b we have 

A N (a,b;T) left = A(a,b;T) + 0(1/N) . (8) 

The problems with the speed of convergence of off- 
diagonal amplitudes in the left prescription can be fixed 
very easily. We find that for all a and b we have 
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e N (V(a)-V(b)) 



A N (a,b-T) left = 

= A(a,b;T) + 0(l/N 2 ) . (9) 



A recent analysis of this method can be found in Jang 
et ali^. Let us note that the crucial step in the deriva- 
tion of the above effective potential from the generalized 
Trotter formula uses the cyclic property of the trace, i.e. 
the above increase in the speed of convergence only holds 
for the partition function and not the amplitudes. A di- 
rect numerical simulation shows that amplitudes calcu- 
lated using this effective potential converge just as fast 
as the amplitudes in the naive left prescription. Said an- 
other way, it is only the integral over all the diagonal 
amplitudes that has the 0(1/N 4 ) behavior and not any 
individual amplitude. A recent investigation by Bond et 
al>i£ has uncovered a 0(1/N 6 ) behavior, however, not for 
the case of a generic theory. At the end let us mention 
that several related investigations dealing with speed of 
convergence have focused on improvements in short time 
propagatio n 16 ! 17 ! 18 or the action 19 . 



III. RELATION BETWEEN DIFFERENT 
DISCRETIZATIONS 



The aim of this and the following section is to present 
a systematic exposition of the relation between different 
discretizations of the same path integral. Throughout 
we will work in the mid-point prescription. We start by 
studying the relation between the 2iV-fold and iV-fold 
discretizations of a given amplitude. From eq. we see 
that we can write the 2./V-fold amplitude as 



A 2N (a, b; T) = 



1 



lTre N 



dqi ■ --dq N -i e 



(12) 
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i.e. in the form of an ./V-fold amplitude given in terms of 
a new action Sn determined by 




d Xl ---dx N e- S2N , (13) 



where S2N is nothing but the 2iV-fold discretization of 
the starting action. In the above formulas we have, for 
convenience, written the 2iV-fold discretized coordinates 
Qo, Qi, . . . , Q2N in terms of q's and x's in the following 
way: Qik — q_k and Q%k-x — %k- Note that we have 
<7o = a, Qn = b, while the N—l remaining q's play the role 
of the dynamical coordinates in the iV-fold discretized 
theory. The x's are the N remaining intermediate points 
that we integrate over in eq. l(T3|) . 

We wish to have Sn belong to the same class of ac- 
tions as Sn- It is not difficult to show that the naively 
discretized action does not satisfy this requirement, i.e. 
the integration of eq. H13(l will yield new types of terms in 
Sn- In fact, the class of actions closed to transformation 



(|13|l is of the form 

Sn = E ( ^ + e NV(q n ) + e N 6l gi {q n ) + 

+ cn &n 92(q n ) + 5 6 n g 3 (q n ) + . . .) . (14) 



The functions appearing in the above effective action also 
depend on the time step ejv- We choose not to display 
this dependence explicitly in order to have a more com- 
pact notation. What is important is that all of these 
functions are regular in the ejv — > limit. Note that 
these effective actions are equivalent to our starting ac- 
tion, i.e. they all have the same continuum amplitudes 
as the starting theory. Using eq. l|13fl and (|14fl one can 
easily derive the following integral relation which deter- 
mines the functions V,gi,g2, ■ ■ ■ in the new action S in 
terms of the related functions in the starting action: 



exp^-e N (v(q n )+6lg 1 (q n )+6t l g2(qn) + .--)^J = (z^) J dy (~ ^-y 2 ^ F (q n + y) , (15) 



where 



-±lnF( X ) = l V (^±^) + lv^ + ^ 



e N v ' 2 V 2 J 2 \ 2 

, 1 1 N2 f gn+1 +x \ 1 2 (x + q r , 

+ 7;Wn+i - x) gi[ + -(x - q n ) g 1 



2 V ^ ' * \ 2 J 2" ' \ 2 

1 

f 



. 1/ ^4 (q n +i+x\ 1 4 (x + q n . 



The above integral equation can be solved for the sim- 
ple cases of a free particle and a harmonic oscillator, and 
gives the well known results. Let us note that the inte- 
gral in eq. I|15|) is in a form that is ideal for an asymptotic 
expansion^, whatever the potential. The time step ejy 
plays the role of a small parameter (in complete parallel 
to the role h plays in the usual semi-classical, or loop, 
expansion of quantum theories). As is usual, the above 
asymptotic expansion is carried through by first Taylor 
expanding F(q + y) around q and then by doing the re- 
maining Gaussian integrals. We find 



(- 

\Tre N 



+00 / 2 

dy exp ( y 2 ) F(q n + y) 

-00 V e N 



E 

771=0 



F (2m) fe ) fe N \ m 
(2m)! 



Finally, from eq. I|15(l and (|17fl we get 

V(q n ) + $1 9\ (<&») + S A n 52(977) + • 



In 



CAT 



E 

ni—0 



(2m)! 



V 4 ) 



(18) 



All that remains is to calculate the F( 2m )(g„)'s using 
eq. Ijltjfl and to expand the potential and all the functions 
gk around the mid-point q n . For this second step we 
make use of the simple relations that follow from the 
definitions of q n and S n - q n +i - q n = q n - q n = <5„/2, 
^77+1 + q n = 2g„ + 5 n /2, and q n + q n = 2q n - S n /2. 
Using these relations we can expand a typical term like 

(Qn+i - qn) g'i ( g " +1 2 +,? " ) to obtain 



( e -f) > ( 17 ) *« > (- 



ySi(«n) + f 9'^) + 



where we have assumed that ejv < 1, i-e. that N > T. 



For example, the expansion of S up to e N is a rather 
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simple exercise. We find the following functional relations 
V = V + e N 



16 y 32 



2048 



128 



9i 



9i 



4 yl 32 



(19) 



:92 



f024 



G4 



92 = ~ 9 2 



6144 



128 



g'l 



Note that in the above relations we expanded V up to 
e N , gi up to ejv an d g 2 up to e N . We also disregarded 
all the higher gk's. The reason for this is that the short 
time propagation of a generic theory satisfies 5„ tx e^r 
while the §k term enters the action multiplied by 5^ k . In 
general, if we expand the new action S to we need to 
evaluate only V (up to e^ 1 ) and the hrst p — 1 functions 
gk (up to e p N ~ 1 ~ k )- Although straight-forward, the task of 
calculating S to large powers of is quite tedious; using 
the symbolic algebra package MATHEMATICA 5.0 we 
have analytically solved the corresponding expressions up 
to p < 9. The memory requirements for this calculation 
grow exponentially with p: the p = 9 calculation used 
just under 2 GB of computer memory. 

At this point it is important to comment on what has 



been achieved so far. Evaluating S to and using this 
new action to calculate the A^-fold discretized amplitude 
An we find 



A N (a, b- T) = A 2N {a, 6; T) + 0(e p N ) 



(20) 



so that, up to 0{e p N ), this amplitude is the same as the 
2A r -fold amplitude calculated with our starting action. 
In this way we have halved the discretization from 2N to 
N. Therefore, a coarser AT- fold discretization using Sn 
does the same job as the 27V-fold discretization of the 
starting theory. In the next section we will consider the 
iteration of this halving procedure. We will derive and 
solve the recursive relations that connect up the 2 s 7V-fold 
and iV-fold discretizations. In particular, we will focus 
on the continuum limit solution when s — » oo, i.e. the 
solution that connects up the continuum theory with its 
iV-fold discretization. 



IV. 



RECURSIVE HALVING 



The iterative process of halving that starts from the 
2 s JV-fold discretization is governed by a recursive rela- 
tion. From the p = 3 case given in eq. 119|) . used in the 
previous section to illustrate the general discretization 
halving scheme, we directly get the sought-after p = 3 
level system of recursive relations 



I4+i 

(ffi)fe+i 
(52W1 



v k 
i 



-k-l 



j(9 



l)k 



32 



:Vl' 



2 2 < 



s-k-l) 



16 



= T(Sl)k + ^Vfc + 



32 



is— k— 1 



:(9 



Y^(92)k 



6144 k 



(4) 



128 (5l) 



2)k 



1024 fe 



(52)* 
1 

64 



32 



1 2 



(91 



2048 k 



(4) 



128 



(9i) k 



(21) 



In the above relations k = 0, 1, 2, . . . , s — 1. The zeroth 
iterate corresponds to the starting action, the last iter- 
ate to the effective action that gives an equivalent iV-fold 
discretization. The £Ar/2 s ~ fe ~ 1 terms represent the time 
step of the fc-th iterate in the discretization halving pro- 
cedure. 



Although the above system of recursive relations is 
non-linear it is in fact quite straight-forward to solve if 
we remember that the system itself was derived via an 
expansion in e^r- Having this in mind we first write all 
the functions as expansions in powers of ejv that are ap- 
propriate to the level p we are working at. In this case 



we are illustrating the procedure for p = 3, so we have 

(ffi)fc = Dk + ^^Ek (22) 
(j92)k = F k ■ 

Putting this into the p = 3 level system of recursive re- 
lations given in eq. (|21|l we find that Ak+i — Af.- Using 
the initial condition Aq = V we find that Ak — V for all 
k. Using this the remaining equations form the following 
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set of linear recursive relations 



2-Bfe+i 
4C fc+1 



C k 





V" 






' 16 






3D>/ E k 


3F fc 


16 


32 2 


4 


V 


' 2 y(4) 




\ 


5 + 512 




V" 






8 

B'l 
8 


n" 

^ + 3F fe 4 
o 


1/(4) 
128 


D'k 






X 4 


384 ' 





W k+1 -D k = — (23) 
r" r>" t/(4) 

8£ fe+ i - £ fc 
16F k+1 - F k 

This system is easily solved for given initial conditions. 
However, what we are really interested in is the contin- 
uum limit solution which is obtained by setting k = s — 1 
in the above expression and taking the limit s — > oo. By 
doing this we are iterating our process of discretization 
halving from the continuum theory down to N. The con- 
tinuum limit solution of the p = 3 level system is simply 

V' 2 vw 

. ~2T 

V" 
~^ 

^(4) 



Vu- 



V + € N 



(51) p=3 

(52) p=3 



V" 

12 + ' 

V-(4) 

480" 



AT 



240 



(24) 



1920 



In the above expressions the label "p — 3" reminds us 
that this is the solution of the continuum limit of the 
recursive relations given in eq. (|21|l describing discretiza- 
tion halving at the p = 3 level. Note that the continuum 
limit solution depends only on the initial potential V, 
i.e. it is not sensitive to initial values of the g k 's as these 
terms all vanish in the continuum limit. In this way we 
have obtained the effective action that gives the best N- 
fold discretization of the starting theory at the p = 3 
level. One can similarly obtain a set of effective actions 
S( p \ one for each value of p. The solution for p = 6 is 
given in the Appendix. Note that each solution contains 
within it all the solutions for lower levels. Solutions for 
larger values of p are a bit more cumbersome, however, 
they are just as easy to use in simulations. Expressions 
up to p — 9 can be found on our web site^i. 

Note that one solves for the continuum limit of the level 
p system of recursive relations but once for all theories, 
i.e. once the solution is found it works for all sufficiently 
smooth potentials V. Actually, the requirement for the 
level p solution is that the starting potential is differen- 
tiable 2p — 2 times. 

The effective action satisfying the continuum limit of 
the discretization halving recursion relations at level p 
leads to an TV- fold amplitude that is equal to the con- 
tinuum amplitude of the starting action up to an 0(e p N ) 
term. Therefore, the continuum limit solution satisfies 



Expectation values can be calculated with the same 
precision using standard discretized estimators, provided 
that the discretized time step on which those observables 
reside is chosen appropriately. For example, the expec- 
tation value of the momentum squared (p 2 (t)) may be 
calculated using the standard estimator 6% if the time 
step from n to n + 1 is shortened to keeping all the 
remaining time steps unchanged. 

The validity of the presented analytical result will be il- 
lustrated in the following section where we present Monte 
Carlo simulations of two different models. To conclude 
this section: we have constructed a general procedure for 
calculating effective actions for any level p. We have 
completed the procedure and found explicit values for the 
effective actions up to and including p = 9. The TV-fold 
amplitudes of the p — 9 effective action differ from the 
continuum amplitudes by a term proportional to 1/N 9 . 



V. NUMERICAL RESULTS AND 
ALGORITHMS 

In this section we illustrate the generic results obtained 
in the previous sections by analyzing the speedup for dif- 
ferent values of p in the case of Monte Carlo simulations 
of two different models. The first model we looked at is 
the anharmonic oscillator with quartic coupling 



V(q) 



1 



A 



(26) 



2 * 4! * 

In Fig. Q] we illustrate how the discretized amplitudes 
Aft tend to the continuum limit for levels p = 1 (naively 
discretized starting action), 2, 4 and 9. The top plot gives 
an overall view of how the discretized amplitudes for a = 
0,6=1 and T = 1 calculated using higher level effective 
actions systematically outperform the ones at lower level. 
In particular we see how they outperform the amplitude 
calculated using the naively discretized action. In the 
bottom plot we show a detail of the top plot which makes 
it easier to qualitatively track the differences between 
amplitudes calculated with effective actions at levels p = 
2,4 and 9. In agreement with eq. 1)25(1 . the curve fitted 
to the p level data is a polynomial in 1/N of the form 



N N p 



C(p) 
NP+ 1 



(27) 



A { v\aMT) = A{aMT) + 0{e p N ) 



(25) 



As derived, all the A^ are (within error bars) equal 
to each other and represent the continuum amplitude 
A(a,b;T) we seek. The continuum value is represented 
in the plots by a dashed line. In all cases the fits were 
done for data with N > 1. The reason that the N = 1 
points were omitted is that for T = 1 we have e\ = 1, i.e. 
these points do not satisfy the condition for asymptotic 
expansion ejv < 1- However, the N — 1 amplitudes are 
quite interesting because they are algebraic expressions 
with no integrations. From the above plot we see that the 
N = 1 amplitudes calculated with higher level effective 
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actions give better and better approximations to the am- 
plitude. We have seen this behavior for other potentials 
as well. 
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FIG. 1: (Color online) (top) Plot of discrete amplitudes ^4^' 
as a function of N for p — 1,2,4 and 9 for an anharmonic 
oscillator with quartic coupling A = 10, time of propagation 
T = 1 from a = to b = 1, N M c = 9.2 • 10 7 . (bottom) Detail 
of the same plot comparing amplitudes for p — 2, 4 and 9. 
In both plots the dashed line represents the continuum limit 
amplitude. 



For T < 1 the N — 1 amplitudes calculated using effec- 
tive actions derived in the previous section represent very 
good algebraic approximations for the case of a general 
theory - larger levels p give better approximations. From 
Fig. ^ we see that these approximations work rather well 
even for the marginal point T = 1. For T > 1 we do 
need to do some integrals numerically in order to get the 
required amplitudes. For high level effective actions it is 
enough to use N = [T] + 1, i.e. to do only [T] integrations 
numerically. 

A quantitative measure of how well the derived effec- 
tive actions perform can be seen in Fig. |2 We see ex- 
plicitly that the p level data differ from the continuum 
amplitudes as polynomials starting with 1/N P . Because 

of this, the deviations from the continuum limit A\ 
become exceedingly small for larger values of p making 




N 



FIG. 2: (Color online) The deviations from the continuum 
limit \Aj^ — A\ as a function of N for p — 1,2,4 and 6 (top to 
bottom). This particular plot is for the case of an anharmonic 
oscillator with quartic coupling A = 10, time of propagation 
T — 1 from a — to b — 1. The number of Monte Carlo 
samples used was Nmc = 9-2 • 10 9 for p — 1,2, Nmc = 
9.2 ■ 10 10 for p = 4, and N M c = 3.68 ■ 10 11 for p = 6. Dashed 
lines correspond to appropriate 1/N polynomial fits to the 
data. The solid lines give the leading 1/N behavior. The 
level p curve has an 1/N P leading behavior. 



it necessary to use ever larger values of Nmc so that 
the Monte Carlo statistical error does not mask these ex- 
tremely small deviations. For p = 6 we see that although 
we used an extremely large number of Monte Carlo sam- 
ples (N M c = 3.68 ■ 10 11 ) the statistical errors become of 
the same order as the deviations already at N > 8. For 
p = 9 this is the case even for N = 2, i.e. we already 
get the continuum limit within a Monte Carlo error of 
around 10~ 8 . 

To make the deviations in Fig.|2visible for large p lev- 
els we needed to run a simulation with very large Nmc- 
This simulation took about a week on our 160 Gflops 
cluster. On the other hand the simulation in Fig. ^ uses 
a much smaller Nmc & n d takes less than one hour to 
complete. In practice we see that the derived effective 
actions give excellent agreement with continuum limit 
amplitudes already for small values of N. Simulations 
with such values of N take a negligible amount of time 
even on a single PC. 

The second model we consider is that of a particle mov- 
ing in a modified Poschl- Teller potential - a well known 
exactly solvable models 



V(q) 



1 a 2 f3((3- 1) 

2 cosh 2 aq 



(28) 



Unlike the anharmonic oscillator the Poschl- Teller poten- 
tial has both a continuous and discrete spectrum. The 
discrete eigenstates have energy 



■(P-l-n) 2 , 



(29) 
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for < n < (3 — 1. Therefore, we see that the model 
has critical values of coupling (3 = 1, 2, 3, . . . at which it 
acquires new bound states. 



0.2644 












0.2642 








p=l 


■ - 


0.264 








P =2 
P =4 


+ 

X 


0.2638 








P =9 




0.2636 












0.2634 












0.2632 




-X - 








- -3ft- - 


- - 




W J 


0.263 


+ 










1 




10 




100 



N 



0.2631600 
0.2631598 h 
0.2631596 
0.2631594 
0.2631592 
0.2631590 







.... *- 




X 










P =2 
P =4 


+ / 

X / 




P =9 


* / 



10 



N 



FIG. 3: (Color online) (top) Plot of discrete amplitudes ^4^' 
as a function of N for p — 1,2,4 and 9 for a particle in a 
modified Poschl- Teller potential with parameters a = 0.5, 
P = 1.5. T = 1, a = 0, b = 1, Nmc = 9.2 ■ 10 7 . (bottom) 
Detail of the same plot comparing amplitudes for p = 2, 4 and 
9. In both plots the dashed line represents the continuum 
limit amplitude. 



Fig. |31 displays how the discretized amplitudes Ap 
tend to the continuum limit for levels p = 1,2,4,9 for 
the potential with a = 0.5 and (3 = 1.5. The same plots 
for the case of a = 0.5 and (3 = 2 (lying on the a critical 
value of (3) are given in Fig. 21 As we can see the effective 
actions work just as well as in the case of the anharmonic 
oscillator. Going through a critical point like (3 = 2 cer- 
tainly affects the physical quantities calculated, however, 
the speedup algorithm is not affected in any way. 

With the increase of p level the complexity of the ex- 
pressions for the effective actions grows exponentially. 
Therefore, the increase in computation time that results 
from using higher p level effective actions also grows ex- 
ponentially as is shown in Fig. [SJ 

As we have seen, by increasing p we drastically improve 




N 



FIG. 4: (Color online) Same plots as in Fig. |5] but for a 
modified Poschl- Teller potential with parameters a = 0.5 and 
(3 = 2. 
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FIG. 5: (Color online) Relative increase in computation time 
that comes about from the increased complexity of expression 
for higher p level effective actions. 



convergence to the continuum limit. An important con- 
sequence of this is that we can obtain the same precission 
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using much smaller values of N, i.e. much coarser dis- 
cretizations. This is at the root of the speedup that we 
find. However, as we have seen, the exponential growth in 
complexity of the effective actions puts an upper bound 
to p levels that can be used. From Fig. [5] we see that 
p = 9 is still far from that upper bound - the gain of 
eight orders of magnitude in the speed of convergence far 
outweighs what is roughly a ten fold increase in compu- 
tation time. 

At the end we briefly comment on two Monte Carlo al- 
gorithms developed for simulations in this section. In the 
first algorithm trajectories are generated by a Gaussian 
distribution function obtained using a semi-classical ex- 
pansion. The computing time of this algorithm scales 
as N 2 ■ Nmc since it is necessary to diagonalize the 
quadratic form in the exponential of the distribution 
function. In the second algorithm we implemented the 
bisection method 5 , which scales as N ■ Nmc- Therefore, 
the bisection algorithm is the method of choice for large 
values of N. On the other hand, our method allows us 
to obtain very precise results using small values of N. 
In that region we have found the two algorithms to be 
comparable both in precision and running time. 

In both algorithms we needed to use a random number 
generator which gives a large number of uncorrelated ran- 
dom numbers in a fashion suitable for parallel program- 
ming. Our primary random number generator was the 
Scalable Parallel Random Number Generator Iibrary22i2i 
(SPRNG). Following the good practise suggested by Fer- 
renberg et al— we have checked all our results using a 
different random number generator. Checks were made 
with the Numerical Recipes' RAN3 generator— with a 
different seed for each MPI process. Agreement was in 
all cases well within a 1-a interval implying that there 
were no hidden systematic errors present in cither the 
algorithms or the random number generators. 

We note in passing that the analytical derivations pre- 
sented in this paper work equally well in both the Eu- 
clidean and Minkowski formalism (with appropriate ie 
regularization), i.e. they are directly applicable to quan- 
tum systems as well as to statistical ones. However, the 
Monte Carlo simulations used to numerically document 
our analytical results necessarily needed to be done in 
the Euclidean formalism. 



VI. CONCLUSION 

To conclude, we have presented an algorithm that leads 
to significant speedup of numerical procedures for the 



calculating of path integrals of a generic theory. The in- 
crease in speed is the result of new analytical input that 
has emerged from a systematic investigation of the re- 
lation between different discretizations of the same the- 
ory. We have presented an explicit procedure for ob- 
taining a set of effective actions that have the same 
continuum limit as the starting action S 1 , but which ap- 
proach that limit ever faster. Amplitudes calculated us- 
ing the A-point discretized effective action sffl satisfy 

A%\a,b]T) = A(a,b;T) + 0(1/N P ), where a and b are 
initial and final states, T the time of propagation, and 
A{a, b; T) the sought-after amplitude of the continuum 
theory. We have obtained and analyzed the effective ac- 
tions for p < 9. In this paper we quote expressions up to 
p = 6 (see the Appendix), the rest can be found on our 
web siteSi. 

At the end we illustrated the obtained generic results 
by analyzing the speedup for different values of p in the 
case of concrete Monte Carlo simulations of two different 
models: anharmonic oscillator with quartic coupling and 
particle in a modified Poschl- Teller potential. 

Extensions of the derived algorithm to M > 1 particles 
and d > I dimensions, as well as to quantum field theo- 
ries are both in progress. In both cases the derivation of 
the analogue of integral eq. l(T3)l does not seem to present 
a problem. The asymptotic expansion used to solve it is 
also directly generalizable. However, the algebraic recur- 
sive relations that determine will be more complex 
and may practically limit us to smaller values of p. 



APPENDIX A: EFFECTIVE ACTION TO p = 6 



In this Appendix we present the effective action at level 
p — 6. Note that this solution contains within it the 
effective actions at all lower levels - all one needs to do 
is to truncate the p = 6 solution at the appropriate order 
in the cat expansion. For example, the effective potential 
V at the p — 3 level is obtained from the p = 6 level 
potential by disregarding term that are and higher. 
Code containing for p < 9 is available on our web 
siteSi. 
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